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ABSTRACT 

Isothermal and adiabatic shocks, which are produced from fast expansion of 
the gas, is simulated with smoothed particle hydrodynamics (SPH). The results 
are compared with the analytic solutions. The algorithm of the program is ex- 
plained and the package, which is written in Fortran, is presented in the appendix 
of this paper. It is possible to change (to complete) the program for a wide variety 
of applications ranging from astrophysics to fluid mechanics. 

Subject headings: Hydrodynamics, methods: numerical, ISM: evolution 

1. Introduction 

Gas dynamical processes are believed to play an important role in the evolution of 
astrophysical systems on all length scales. Smoothed particle hydrodynamics (SPH) is a 
powerful gridless particle method to solve complex fluid-dynamical problems. SPH has a 
number of attractive features such as its low numerical diffusion in comparison to grid based 
methods. An adequate scenario for SPH application is the unbound astrophysical problems, 
especially on the shock propagation (see, e.g., Liu & Liu 2003). In this way, the basic 
principles of the SPH is written in this paper and the simulation of isothermal and adiabatic 
shocks are applied to test the ability of this numerical simulation to produce known analytic 
solutions. 

The program is written in Fortran and is highly portable. This package supports only 
calculations for isothermal and adiabatic shock waves. It is possible to change (to complete) 
the program for a wide variety of applications ranging from astrophysics to fluid mechanics. 
The program is written in modular form, in the hope that it will provide a useful tool. I ask 
only that: 

• If you publish results obtained using some parts of this program, please consider ac- 
knowledging the source of the package. 
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• If you discover any errors in the program or documentation, please promptly commu- 
nicate the to the author. 



2. Formulation of Shock Waves 

An extremely important problem is the behavior of gases subjected to compression 
waves. This happens very often in the cases of astrophysical interests. For example, a small 
region of gas suddenly heated by the liberation of energy will expand into its surroundings. 
The surroundings will be pushed and compressed. Conservation of mass, momentum, and 
energy across a shock front is given by the Rankine-Hugoniot conditions (Dyson & Williams 
1997) 

pivi = p 2 v 2 (1) 
Pl vl + Kp\ = p 2 v 2 2 + Kp\ (2) 

\v\ + ^-rKpJ- 1 = \v\ + ^-Kpl 1 + Q (3) 
Z 7 — 1 Z 7 — 1 

where the equation of state, p = Kp 1 , is used. In adiabatic case, we have Q = 0, and for 
isothermal shocks, we will set 7 = 1. 

We would interested to consider the collision of two gas sheets with velocities v in the 
rest frame of the laboratory. In this reference frame, the post-shock will be at rest and the 
pre-shock velocity is given by v± = v + 1> 2 , where t> 2 is the shock front velocity. Combining 
equations (l)-(3), we have 



^2 = ao["2 + V 1 + J + (7 - - ?)] (4) 



where ao = lKp[ 1 is the sound speed, M = i>o/ao is the Mach number, b and q are defined 

as 

Substituting (4) into equation (1), density of the post-shock is given by 



(6) 
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3. SPH Equations 

The smoothed particle hdrodynamics was invented to simulate nonaxisymmetric phe- 
nomena in astrophysics (Lucy 1977, Gingold & Monaghan 1977). In this method, fluid is 
represented by N discrete but extended/smoothed particles (i.e. Lagrangian sample points). 
The particles are overlapping, so that all the physical quantities involved can be treated 
as continues functions both in space and time. Overlapping is represented by the kernel 
function, W a b = W(r a — r b ,h a b), where h a b = (h a + hb)/2 is the mean smoothing length 
of two particles a and b. The continuity, momentum and energy equation of particle a 
are (Monaghan 1992) 

Pa = m b W ab (7) 

b 

^ = - E ^ + - + Kab)V a W ab (8) 
at ^ pa Pb 

at 2 ^ p a p b 

where v a6 = v a - v 6 and 

u ab = { — — ' 11 Vab - Tab < U ' (10) 

1 0, otherwise, 

is the artificial viscosity between particles a and b, where p ab = \{p a + Pb) is an average 
density, a ~ 1 and (3 ~ 2 are the artificial coefficients, and /i a & is defined as its usual form 

with ~ 0. 1 and fo a f, = \{h a + The signal velocity, v S i g , is 

^i 9 = ^(Ca + C 6 ) (12) 

where c a and c 6 are the sound speed of particles. The SPH equations are integrated using 
the smallest time-steps 

At a = C cour MIN[^- (r^r, r-T^n, TtAttT' TT%|] ( 13 ) 
| v I I a_i | | du a /dt | | dh a /dt \ | dp a /dt \ 

where C cour ~ 0.25 is the Cour ant number. 
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4. Results and Prospects 

The chosen physical scales for length and time are [/] = 3.0 x 10 16 m and [t] = 3.0 x 10 13 s, 
respectively, so the velocity unit is approximately lkm. s -1 . The gravity constant is set 
G — 1 for which the calculated mass unit is [m] = 4.5 x 10 32 kg. There is considered two 
equal one dimensional sheets with extension x = 0. 1 [Z] , which have initial uniform density 
and temperature of ~ 4.5 x 10 8 m~ 3 and ~ 10K, respectively. 

Particles with a positive x-coordinate are given an initial negative velocity of Mach 5, 
and those with a negative x-coordinate are given a Mach 5 velocity in the opposite direction. 
In adiabatic shock, with M = 5, the post-shock density must be 2.9, which is obtained from 
analytic solution (6) with Q = and 7 = 2. The Results of adiabatic shock are shown in 
Fig. 1-3. In isothermal shock, with M = 5, the post-shock density must be 26.9, which is 
obtained from analytic solution Equ. (6) with 7 = 1. The Results of isothermal shock are 
shown in Fig. 4-5. Algorithm of the program is shown in Fig. 6. 
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Fig. 1. — The density of adiabatic shock, with M = 5, Q = 0, and 7 = 2. 
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Fig. 2. — The velocity of adiabatic shock, with Mq = 5, Q = 0, and 7 = 2. 
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Fig. 3. — The temperature of adiabatic shock, with M = 5, Q = 0, and 7 = 2. 
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Fig. 4. — The density of isothermal shock, with Mq = 5 and 7 = 1. 
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Fig. 5. — The velocity of isothermal shock, with Mq = 5 and 7 = 1. 
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SHOCKSET 

set initial values of the adiabatic or 
isothermal shock waves 






SHOCKEND 

check the end conditions of the adiabatic 
or isothermal shock waves 


1 — > 


I 




TREENNSL 

make tree, find nearest neighbors, and the 
smoothing lengths 


i 




STATES 

find density, pressure, temperature, energy, 
sound speed, and different rates 






ADVANCE 

find minimum time- step and advance 
the particles at one time-step 





Fig. 6. — Algorithm of the smoothed particle hydrodynamics for simulation of isothermal 
and adiabatic shocks. 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
This program is provided to simulate the adiabatic 
and isothermal shock waves 

********************** 

Mohsen Ne jad-Asghar 
nasghar@dubs . ac . ir 
January 2006 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
PROGRAM SHOCK 
INCLUDE 1 par am . inc ' 

PRINT*, 'isothermal or adiabatic shock?' 

PRINT*, 1 adiabatic=l 1 

PRINT*, 1 isothermal-2 ' 

READ ( * , * ) isorad 

CALL SHOCKSET 

CALL SCENARIOS 

! investigate the end condition of simulation 
10 CALL SHOCKEND 

! advance system at one time-step 
CALL ADVANCE 
GOTO 10 

END PROGRAM SHOCK 



//////////////////////////////////////////////////////////////////// 



SUBROUTINE SHOCKSET 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
This subroutine generates initial particle information for 
adiabatic or isothermal one dimensional shock 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE ' par am. inc 1 
REAL extx, delx, mach 
nbody-4 00 
IF (nbody > maxn) 
& CALL TERROR (' SHOCKSET : SPH number is very large') 

mxcell=nsubc* nbody 
node=nbody+mxcell 
incell=nbody+l 

! units of length, time, and mass 
ul=3 . 0el6 
ut=3 . 0el3 
um=4 . 5e32 

! extension of each sheet in x direction (ul) 
extx=0 . 1 

! positions of SPH particles 
delx=2 . 0*extx/nbody 
DO i=l, nbody 

pos (i, 1) =extx-i*delx+delx/2 . 
END DO 

! density of SPH particles (um/ul A 3) 
den=l . 

! temperature of SPH particles (K) 
temp=10 . 

! masses of SPH particles (um) 
DO i=l, nbody 

mass (i) =delx*den (i) 
END DO 

! molecular weight relative to the mass of hydrogen 
xmu=2 . 

! hydrogen mass 
xmh=l . 67e-27 

! Bolt zman constant 
xkb=l . 38e-23 

xKK= (xkb/ (xmu*xmh) ) / (ul/ut) **2 
IF (isorad == 1) THEN 

! polytropic index (adiabatic case) 

gamma=2 . 

! energy of SPH particles 

u=xKK*temp/ (gamma-1) 
ELSEIF {isorad ==2) THEN 

! polytropic index {isothermal case) 

gamma=l . 
END IF 

! sound speed 

sound=SQRT (gamma*xKK*temp) 
! Mach number 
mach=5 . 

IF (isorad == 1) THEN 

! relative density after simulation for adiabatic case 
cons=SQRT (16 . 0+ ( {3-gamma) **2+8* (gamma-1) ) *mach*mach) 
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denfinal= (const ( gamma+1 ) *mach) / (cons- ( 3 -gamma) *mach) 
ELSEIF (isorad == 2 ) THEN 

! relative density after simulation for isothermal case 

cons=SQRT ( 4 . 0+mach*mach) 

denf inal- (cons+mach) / (cons-mach) 
END IF 

PRINT*, 'relative density after simulation must be', denfinal 
! velocity of SPH particles 
DO i=l, nbody 

IF (pos (i, 1) > 0.0)THEN 
vel (i, 1 ) =- sound (i ) *mach 

ELSE 

vel (i, 1 ) =+ sound (i ) *mach 
ENDIF 
END DO 

! initial smoothing lengths 
diminv=l . 0/dim 
DO i=l, nbody 

hh (i) =2 . 0* (mass (i) /den (i) ) **diminv 
END DO 
tnow-0 . 

PRINT*, 'setup is successfully completed' 
PAUSE 'press ENTER to continue' 
END SUBROUTINE SHOCKSET 
! ==================================================================== 

\ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
!==================================================================== 

SUBROUTINE SHOCKEND 
Icccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
! This subroutine find the end conditions for adiabatic one 

! dimensional shock 

I cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE 'param.inc' 
INTEGER nrev 
nrev=0 

DO i=l, nbody 12 
IF(vel(i,l) > 0.0) nrev=nrev+l 

DO i=nbody/2+l, nbody 

IF(vel(i,l) < 0.0) nrev=nrev+l 
END DO 

PRINT*, 'tnow=', tnow, ' max density=', MAXVAL (den) 

! stop the program when the reflection waves occur 
IF (nrev > nbody/ 10 ) THEN 
CALL SAVEFIG 

PRINT*, 'reversed particles=', nrev 
PAUSE 
ENDIF 

END SUBROUTINE SHOCKEND 
!==================================================================== 

\ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
!==================================================================== 

SUBROUTINE SAVEFIG 
Icccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
! This program writes particle information into different 

! files to draw the figures 

Icccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE 'param.inc' 
OPEN (1, file=' posden.dat ' ) 
OPEN (2, file=' pospre.dat ' ) 
OPEN (3, file=' posvel.dat ' ) 
OPEN(4, file='posu.dat' ) 
DO i=l, nbody 

WRITE (1,*) pos(i,l), den(i), tnow 

WRITE (2,*) pos(i,l), pre(i), tnow 

WRITE (3,*) pos(i,l), vel(i,l), tnow 

WRITE(4,*) pos(i,l), temp(i), u(i), tnow 
END DO 
CLOSE (1) 
CLOSE (2) 
CLOSE (3) 
CLOSE (4) 
END SUBROUTINE SAVEFIG 



\ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
!==================================================================== 

SUBROUTINE SCENARIOS 
I cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
! This subroutine switches for different scenarios in simulation 

Icccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
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INCLUDE 'param.inc' 

! skf--> smoothing kernel function? 

! =1 : Gauss kernel (Gingold & Monaghan 1981) 

! =2: spline-base kernel (Monaghan 1985) 

! =3: Quintic kernel (Morris 1997) 

skf=2 

! nnssl — > nearest neighbors and smoothing length? 

! =1 : fixed smoothing length 

! =2 : variable smoothing length 

nnssl=l 

! dsm--> density summation method? 

! -1: summation model without continuity 

! =2: use continuity equation 

dsm=l 

! the artificial shear viscosity? 
alphas=l . 

! the artificial bulk viscosity? 
betas=2 . 
END SUBROUTINE SCENARIOS 



1 1 i 1 1 1 1 1 i 1 1 i 1 1 i 1 1 1 1 1 i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 i 1 1 i 1 1 i 1 1 1 1 1 i 1 1 i 1 1 i 1 1 1 1 1 i 1 1 i 1 1 i 



SUBROUTINE ADVANCE 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

This subroutine advances the particles at one time-step 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE 'param.inc' 
REAL velO (nbody, dim) 

REAL denO (nbody) , hhO (nbody) , uO (nbody) 
! advance particles at first time-step 
IF (tnow == 0.0) THEN 

! make tree and find neighbors, smoothing 

! length, and density 

CALL TREENNSL 

! find all states of the system 
CALL STATE 

! find minimum time-step 

CALL COURANT 

DO 1=1, nbody 

IF (nnssl == 2) hh (i) =hh (i) thhdot (i) *dtmin/2 . 
IF(dsm == 2) den (i) =den (i) tdendot (i) *dtmin/2 . 
IF(isorad == 1) u ( i ) =u ( i ) tudot ( 1 ) *dtmin/2 . 
DO j = l, dim 

vel (i, j ) =vel (i,j)+acc(i,j) *dtmin/2 . 
pos (i, j ) =pos (i, j) +vel (i, j) *dtmin 
END DO 

END DO 

tnow=tnow+dtmin 
RETURN 
END IF 

! advance particles at first half time-step 
DO i=l, nbody 
hhO (i) =hh (i) 

IF (nnssl == 2) hh (i) =hh (i) thhdot (1) *dtmin/2 . 
denO { i ) =den ( 1 ) 

IF(dsm == 2) den (i) =den (1) tdendot (i) *dtmin/2 . 
uO (i)=u(i) 

IF(isorad == 1) u ( 1 ) =u ( i ) tudot ( i ) *dtmin/2 . 
DO j=l, dim 

velO (i, j) =vel (1, j) 

vel (i, j ) =vel (1, j) tacc (i, j) *dtmin/2 . 
END DO 
END DO 

dtminl=dtmin 

! make tree and find neighbors, smoothing 
! length, and density 
CALL TREENNSL 

! find all states of the system 
CALL STATE 

! find minimum time-step 
CALL COURANT 

dtmin2=dtminl/2 . Otdtmin/2 . 

! advance particles at second half time-step 
DO i=l, nbody 

IF (nnssl == 2) hh ( 1 ) =hh0 ( i ) thhdot ( i ) *dtmin2 
IF(dsm == 2) den (i) =den0 (i) tdendot (i) *dtmin2 
IF(isorad == 1) u (i) =u0 (i) tudot (i) *dtmin2 
DO j=l, dim 

vel (i, j) =vel0 (i, j) tacc (i, j) *dtmin2 
pos (i, j)=pos (1, j) tvel (1, j) Mtmin 
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END DO 
END DO 

tnow=tnow+dtmin 

END SUBROUTINE ADVANCE 



1 1 i 1 1 1 1 1 i 1 1 i 1 1 i 1 1 1 1 1 i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 i 1 1 i 1 1 i 1 1 1 1 1 i 1 1 i 1 1 i 1 1 1 1 1 i 1 1 i 1 1 i 



SUBROUTINE COURANT 
Icccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
! This subroutine evaluates time-step for each particle 

Icccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE 'param.inc' 

REAL dtl, dt2, dt3, dt4, dt5, delt(5) 
REAL accO, velO 
REAL dt (nbody) 
DO i=l, nbody 

velO=0 . 

accO=0 . 

DO k=l, dim 

velO=velO+vel (i,k) **2 
acc0=acc0+acc(i,k)**2 

END DO 

velO=SQRT (velO) 

accO=SQRT (accO) 

dtl=0 . 

dt2=0 . 

dt3=0 . 

dt4=0 . 

dt5=0 . 

IF(velO /= 0.0) dtl=hh (i) /velO 

IF(accO /= 0.0) dt2=SQRT (hh (i) /accO) 

IF(nnssl == 2 .AND. hhdot(i) /= 0.0) 
s dt 3=hh ( i ) / ABS {hhdot { i ) ) 

IF(dsm == 2 .AND. dendot(i) /= 0.0) 
& dt4=den (i) /ABS (dendot (i) ) 

IF(udot(i) /= 0.0) dt5=u (i) /ABS (udot (i) ) 

j0=0 

IF (dtl /= 0.0) THEN 

j0= j0+l 

delt ( j0)=dtl 
ENDIF 

IF(dt2 /= 0.0) THEN 

j0= j0+l 

delt ( j0)=dt2 
ENDIF 

IF(dt3 /= 0.0)THEN 

j0= j0+l 

delt ( j0)=dt3 
ENDIF 

IF(dt4 /= 0.0) THEN 

j0= j0+l 
delt ( j0)=dt4 
ENDIF 

IF(dt5 /= 0.0) THEN 

j0= j0+l 

delt ( j0)=dt5 
ENDIF 

dt ( i ) =MAXVAL ( delt ) 
DO j=l, jO 

dt (i) =MIN (dt (i) , delt ( j) ) 
END DO 
END DO 

dtmin=cour*MINVAL (dt) 

IF(dtmin == 0.0) CALL TERROR (' COURANT : zero time-step') 
END SUBROUTINE COURANT 



1 1 i 1 1 1 1 1 i 1 1 i 1 1 i 1 1 i 1 1 i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 i 1 1 i 1 1 i 1 1 1 1 1 i 1 1 i 1 1 1 1 1 1 1 1 i 1 1 i 1 1 / 



SUBROUTINE TREENNSL 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

This subroutine makes tree, finds sorted nearest neighbors, 

and estimates appropriate smoothing length and density of 

all particles self-consistently . 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE 'param.inc' 

! construction of tree according to J.E. Barnes 
CALL TREE 

! find nearest neighbors and smoothing lengths 
CALL NNSL 
END SUBROUTINE TREENNSL 
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! ==================================================================== 

\ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
!==================================================================== 

SUBROUTINE TREE 

Icccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
! This subroutine constructs tree, finds its properties, and 

! evaluates the gravitational acceleration 

Icccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE 'param.inc' 
! check particle overlapping 
DO i=l, nbody-1 
DO j=i+l,nbody 
ri j=0 . 
DO k=l, dim 

rij=rxj+(pos(i,k)-pos(j,k) ) * (pos ( i , k) -pos ( j , k) ) 
END DO 

IF(rij == 0.0) CALL TERROR (' TREE : particle overlapping') 
END DO 
END DO 

! construct the octal-tree body-by-body 
CALL TREELOAD 

! find the properties of bodies and cells 
CALL TREEPROP 
END SUBROUTINE TREE 
!==================================================================== 

\ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
!==================================================================== 

SUBROUTINE TREELOAD 
Icccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
! This subroutine constructs the octal-tree body-by-body 

Icccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

INCLUDE 'param.inc' 

INTEGER MKCELL, p 

! allocate root cell 

ncell=0 

root=MKCELL ( ) 

! expand size of the root cell to hold all bodies 
dist=2 . 05*MAXVAL (ABS (pos) ) 
clsize (root) =dist 

! store geometric midpoint of root cell 
DO i =1, dim 

mid(root,i)=0.0 
END DO 

! load bodies into the new tree, one at a time 
DO p=l, nbody 

CALL LDBODY(p) 
END DO 

END SUBROUTINE TREELOAD 

\ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 / 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



SUBROUTINE LDBODY(p) 
Icccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
! This subroutine loads body p into tree structure 

Icccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE 'param.inc' 

INTEGER p, q, qind, SBINDX, MKCELL, c, pO 

! set {q, qind) pair in correct subcell of root 

q=root 

qind=SBINDX (p, q) 

! loop descending tree until an empty subcell is found 
10 IF (subp (q, qind) /= ) THEN 

! if subp {q, qind) is a body, extend the tree with a new cell 
IF (subp (q, qind) < incelDTHEN 

! allocate an empty cell to hold both bodies 
c— MKCELL ( ) 

! locate midpoint of new cell 
DO i=l, dim 

IF (pos (p, i) >= mid (q, i) ) THEN 

mid(c, i) =mid (q, i) tclsize (q) /4 . 
ELSE 

mid(c, i) =mid (q, i) -clsize (q) /4 . 
END IF 
END DO 

! set size of new cell 
clsize (c) =clsize (q) 12 . 

! store old body in appropriate subcell within new cell 
p0=subp (q, qind) 
subp (c, SBINDX (pO, c) ) =p0 
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! link new cell into tree in place of old body 
subp (q, qind) =c 
ENDIF 

! at this point, the node indexed by (q, qind) is known to 
! be a cell, so advance to the next level of tree, and loop 

q=subp (q, qind) 

qind-SBINDX (p, q) 

GOTO 10 
ENDIF 

! found place in tree for p, so store it there 
subp (q, qind) =p 
END SUBROUTINE LDBODY 



//////////////////////////////////////////////////////////////////// 



INTEGER FUNCTION MKCELL { ) 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

This function allocates a cell and returns its index 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE 'param. inc 1 

! check remaining space for a new cell 
IF(ncell > mxcell) THEN 

CALL TERROR ( 'MKCELL: no more memory') 
ENDIF 

! increment cell counter, initialize new cell pointer 

ncell=ncell+l 

MKCELL=ncell+nbody 

! zero pointers to subcells of new cell 
DO i=l , nsubc 

subp {MKCELL, i) =0 
END DO 
END FUNCTION MKCELL 



! //////////////////////////////////////////////////////////// //////// 



INTEGER FUNCTION SBINDX{p,q) 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

This function computes subcell index for node p within cell q 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

INCLUDE ' param. inc 1 

INTEGER p, q 

! initialize subindex to point to lower left subcell 
SBINDX=1 

! loop over all spatial dimensions 
DO i=l, dim 

IF (pos (p, i) >= mid (q, i) ) SBINDX=SBINDX+2 * * (dim-i) 
END DO 
END FUNCTION SBINDX 



//////////////////////////////////////////////////////////////////// 



SUBROUTINE TREEPROP 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
This subroutine checks tree structure, assigns critical radius 
for each cell, computes cell masses , cm. positions, and 
quadrupole moments 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE ' param. inc 1 
INTEGER p, q 
REAL posO(dim), dist2 

! list cells in order of descending size 
CALL SORTLIST 

! loop processing cells from smallest to root 
DO i=ncell, 1, -1 
p=sortind ( i) 

! check that p is a cell 

IF (p < incell) CALL TERROR (' TREEPROP : wrong cell') 
! zero accumulators for this cell 
mass (p) =0 . 
pos0=0 . 

! compute cell properties as sum of properties 
! of its subcells 
DO j=l, nsubc 
q=subp (p, j) 

! only access cells which exist 
IF (q /= 0) THEN 

! sum properties of subcells to obtain 

! values for cell p 

mass (p) =mass (p) +mass (q) 

DO k=l, dim 



posO (k)=posO (k) +mass (q) *pos (q, k) 
END DO 
ENDIF 
END DO 

! normalize center of mass coordinates by total cell mass 
DO j=l, dim 

posO (j) =posO (j) /mass (p) 
END DO 

! check tree, compute cm-to-mid distance 
! and assign cell position 
dist2=0 . 
DO j=l, dim 

IF(posO(j) < mid (p, j) -clsize (p) 12 . .OR. 
& posO(j) >= mid (p, j ) +clsize (p) 12 . 0) 

s CALL TERROR (' TREEPROP : tree structure error') 

dist2=dist2+ (posO ( j) -mid(p, j) ) **2 

! copy cm position to cell. This overwrites the midpoint 
pos (p, j) =pos0 ( j ) 
END DO 

! assign critical radius for cell, adding offset 
! from midpoint for more accurate forces. This 
! overwrites the cell size 

rcrit2 (p) = (clsize (p) /theta+SQRT (dist2) ) **2 
! compute quadrupole moments 
DO j=l, nquad 

quad (p, j) =0 . 
END DO 

! loop over descendants of cell p 
DO j=l, nsubc 
q=subp (p, j) 
IF (q /= 0) THEN 

! sum properties of subcell q to 
! obtain values for cell p 
DO m=l, MIN (2, dim) 
DO n=m, dim 

1= (m-1) * (dim-1) +n 

quad (p, 1) =quad (p, 1) +3 . 0*mass (q) * (pos (q,m) 
& -pos (p, m) ) * (pos (q, n) -pos (p, n) ) 

IF (m == n) THEN 
DO k=l, dim 

quad (p, 1) =quad(p, 1) -mass (q) * (pos (q,k) 
& -pos (p, k) ) **2 

END DO 
ENDIF 

! if q itself is a cell, add its moments too 
IF (q >= incell) quad (p, 1 ) =quad (p, 1 ) tquad (q, 1 ) 
END DO 
END DO 
ENDIF 
END DO 
END DO 

END SUBROUTINE TREEPROP 
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SUBROUTINE SORTLIST 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

This subroutine sorts cells from largest (root) to smallest 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

INCLUDE 'param.inc' 

INTEGER facell, lacell, nacell 

! start scan with root as only active cell 

sortind (1 ) =root 

facell=l 

lacell=l 

! loop while active cells to process 
10 IF (facell <= lacell) THEN 

! start counting active cells in next iteration 
nacell=lacell 

! loop over subcells of each active cell 
DO i=l, nsubc 

DO j=facell, lacell 

! add all cells on next level to active list 
IF (subp (sortind ( j) , i) >= incell ) THEN 
nacell=nacell+l 
IF (nacell > maxn*nsubc) 
& CALL TERROR (' SORTLIST : overflow') 

sortind (nacell) =subp (sortind(j),i) 
ENDIF 
END DO 
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END DO 

! advance first and last active cell indices, and loop 
f acell=lacell+l 
lacell=nacell 
GOTO 10 
END IF 

! above loop should list all cells; check the count 
IF(nacell /= ncell) THEN 
WRITE (*,*) nacell, ncell 

CALL TERROR (' SORTLIST : inconsistent cell count') 
END IF 

END SUBROUTINE SORTLIST 
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SUBROUTINE NNSL 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
This subroutine finds the nearest neighbors and smoothing 
length of particles 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE 'param.inc' 
REAL DIST 
INTEGER p, q, qsub 
IF(nnssl == 1 ) THEN 

! fixed smoothing length 
DO p=l, nbody 

hh (p) =1.5* (mass (p) /den (p) ) ** (1.0/dim) 
! effective radius according to hh 
IF(skf == 1) rp=3.0*hh(p) 
IF(skf == 2) rp=2.0*hh(p) 
IF(skf == 3) rp=3.0*hh(p) 
neighb (p) =0 

! loop processing cells from root to smallest 
DO i=l, ncell 
q=sort ind { i ) 
rr=rp+clsize (q) 
IF (DIST (p, q, rr, 0) < 0.0)THEN 
! accepted: permit descent 
DO j=l, nsubc 
qsub=subp (q, j ) 

IF (qsub /= .AND. qsub < incell) THEN 
! a body: skip self-consideration 
IF (qsub /= p) THEN 
! test its spacing 
rr=rp 

IF (DIST (p, qsub, rr, 1) < 0.0) THEN 

! accepted as a nearest neighbor 
neighb (p) ^neighb (p) +1 
! check number of nearest neighbors 
IF (neighb (p) == nbody) 
S CALL TERROR ( 'NNSL: too many') 

neighblist (p, neighb (p) ) =qsub 
ENDIF 
END IF 
ENDIF 
END DO 
ENDIF 
END DO 

CALL SORTNEIGHB (p) 
END DO 
ELSEIF(nnssl == 2 ) THEN 

! variable smoothing length 
DO p=l, nbody 
numiter=0 
10 numiter=numiter+l 

IF(numiter > 20) THEN 

WRITE!*,*) p, dennew, hhnew 
pause 

CALL TERROR) 'NNSL: too many iteration') 
ENDIF 

! first use the smoothing length at this time, which 

! is advanced via dh/dt= (h/dim) *divvel and find 

! effective radius according to this hh 

IF(skf == 1) rp=3.0*hh(p) 

IF(skf == 2) rp=2.0*hh(p) 

IF(skf == 3) rp=3.0*hh(p) 

neighb (p) =0 

! loop processing cells from root to smallest 
DO i=l, ncell 
q=sortind (i) 
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rr=rp+clsize (q) 

IF (DIST (p, q, rr, 0) < 0.0)THEN 

! accepted: permit descent 

DO j=l, nsubc 
qsub-subp (q, j ) 

IF(qsub /= .AND. qsub < incell) THEN 
! a body: skip self-consideration 
IF (qsub /= p) THEN 
! test its spacing 
rr=rp 

IF (DIST (p, qsub, rr, 1) < 0.0) THEN 

! accepted as a nearest neighbor 
neighb (p) -neighb (p) +1 
! check number of nearest neighbors 
IF {neighb (p) == nbody) 
S CALL TERROR ( 'NNSL: too many') 

neighblist (p, neighb (p) ) =qsub 
ENDIF 
END IF 
ENDIF 
END DO 
ENDIF 
END DO 

! next find density by a summation over the particles 
hmin-hh (p) 

dennew-mass (p) *W (p, p) 
DO jcursor=l, neighb (p) 

j=neighblist {p, jcursor) 

dennew-dennew+mass { j ) *W (p, j ) 
END DO 

! change the smoothing length via the 
! proportionally { 1/den) A ( 1 /dim) 
hhnew-2 . * (mass (p) / dennew) * * ( 1 . /dim) 
! check convergence of smoothing length 
hf rac=ABS (hhnew-hh (p) ) /hh (p) 
IF(hfrac > 0.01) THEN 
hh (p) =hhnew 
GOTO 10 
ENDIF 

CALL SORTNEIGHB (p) 
END DO 
ENDIF 
END SUBROUTINE NNSL 
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FUNCTION DIST (i, q, rr, mode) 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
This function estimates the spacing criterion between particle 
p and node q 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE 'param.inc' 
REAL DIST, rpq 
INTEGER q 
rpq=0 

IF (mode == 1 ) THEN 
DO j=l, dim 

a=pos (q, j ) -pos (i, j ) 

rpq=rpq+a*a 
END DO 

DIST=rpq-rr*rr 
ELSE 

DO j=l, dim 

a=ABS (pos (q, j) -pos (i, j) ) 
rpq-MAX ( rpq , a ) 
END DO 
DIST=rpq-rr 
ENDIF 
END FUNCTION DIST 
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SUBROUTINE SORTNEIGHB ( i ) 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
This subroutine sorts the nearest neighbors at ascending 
distance to particle i 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE 'param.inc' 
REAL rij (neighb (i) ) 

INTEGER indx (neighb (i) ) , indxn (neighb (i) ) 
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ri j=0 . 

DO jcursor-1, neighb(i) 
j=neighblist (i, jcursor) 
DO k=l, dim 

rij(jcursor)=rij(jcursor) + (pos(i,k)-pos(j,k))* A 2 
END DO 

rij (jcursor) =SQRT (rij (jcursor) ) 
END DO 

CALL INDEXX (neighb (i) , ri j, indx) 
DO j=l, neighb (i) 

indxn ( j ) =neighblist { i , indx ( j ) ) 
END DO 

DO j=l, neighb (i) 

neighblist (i, j ) =indxn ( j ) 
END DO 

END SUBROUTINE SORTNEIGHB 



SUBROUTINE INDEXX (n, arr, indx) 
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
! This subroutine indexes an arry arr(l:n), i.e. output the 

! array indx(l:n) such that arr(indx(j)) is in ascending order 

! for j=l,2,..,n. According to 'Numerical Recipes', Press et al. 

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

INCLUDE 'param.inc' 

INTEGER n, indx (n) ,M, NSTACK 

REAL arr (n) 

PARAMETER (M=7 , NSTACK=50 ) 

INTEGER i, indxt, ir, itemp, j, jstack, k, 1, istack (NSTACK) 
REAL a 
DO j=l,n 

indx ( j) = j 
END DO 
jstack=0 
1=1 
ir=n 

1 IF (ir-1. It .M) THEN 

DO j=l+l,ir 
indxt=indx ( j ) 
a=arr (indxt) 
DO i=j-l,l,-l 

IF (arr (indx (i) ) .le.a) GOTO 2 

indx (i+1) =indx (i) 
END DO 
i=l-l 

2 indx (i+1) =indxt 
END DO 

IF (jstack . eq. 0) RETURN 

ir=istack (jstack) 
l=istack(jstack-l) 
jstack=jstack-2 
ELSE 

k= (1+ir) 12 
itemp=indx (k) 
indx (k) =indx (1+1) 
indx ( 1+1 ) =itemp 

IF (arr (indx (1) ) . gt . arr (indx (ir) ) ) THEN 

itemp=indx ( 1 ) 

indx (1) =indx (ir) 

indx ( ir ) =itemp 
ENDIF 

IF (arr (indx (1+1) ) . gt . arr (indx (ir) ) ) THEN 

itemp=indx ( 1+1 ) 

indx ( 1+1 ) =indx ( ir ) 

indx (ir) =itemp 
ENDIF 

IF (arr (indx (1) ) . gt . arr (indx (1+1) ) ) THEN 

itemp=indx ( 1 ) 

indx (1) =indx (1+1) 

indx ( 1+1 ) =itemp 
ENDIF 
1=1+1 
j=ir 

indxt =indx ( 1 + 1 ) 
a=arr ( indxt ) 

3 CONTINUE 
l=i+l 

IF (arr (indx (1) ) .It. a) GOTO 3 
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j=j-i 

IF (arr (indx ( j) ) .gt.a) GOTO 4 
IF(j.lt.i) GOTO 5 
itemp=indx ( i ) 
indx { i ) =indx ( j ) 
indx ( j ) =itemp 
GOTO 3 

5 indx (1+1) =indx ( j) 

indx { j ) =indxt 
j stack- jstack+2 

IF ( jstack . gt . NSTACK) PAUSE ' NSTACK too small in indexx ' 
IF (ir-i+1 . ge . j-1) THEN 

istack (jstack) =ir 

istack(jstack-l)=i 

ir= j-1 
ELSE 

istack (jstack) =j-l 
istack ( jstack-1) =1 
l = i 
END IF 
END IF 
GOTO 1 
END SUBROUTINE INDEXX 



1 1 i 1 1 1 1 1 i 1 1 i 1 1 i 1 1 1 1 1 i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 i 1 1 1 1 1 i 1 1 1 1 1 i 1 1 1 1 1 1 1 1 1 1 1 i 1 1 i 1 1 i 



SUBROUTINE STATE 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

This subroutine finds different states 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

INCLUDE 'param.inc' 

IF(dsm == 1) CALL DENSUM 

CALL DIVVELO 

! find density rate and smoothing length rate 
IF {dsm == 2) dendot--den A divvel 
! find smoothing length rate 
IF(nnssl == 2) hhdot= (hh/dim) *divvel 

! find pressure, sound speed, and energy of particles 
CALL PRESOUNDENG 

! find time-rate of velocity (acceleration) and 
! energy of particles 
CALL RATES 
END SUBROUTINE STATE 



1 1 1 1 1 1 1 1 1 1 1 i 1 1 1 1 1 i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 i 1 1 i 1 1 i 1 1 1 1 1 1 1 1 i 1 1 1 1 1 1 1 1 1 1 1 i 1 1 i 



SUBROUTINE DENSUM 

I cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
! This subroutine estimates the density via normalization 

! summation method 

Icccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE 'param.inc' 
REAL sumW(nbody) 

! firstly integration of the kernel 
! secondly density integration 
DO i=l, nbody 
hmin=hh ( i ) 

den ( i) =mass (i) *W (i, i) 
DO jcursor=l, neighb(i) 

j=neighblist (i, jcursor) 

hmin= (hh (i) +hh ( j) ) 12 . 

den (i) =den ( i) +mass ( j ) *W (i, j ) 
END DO 
END DO 

DO i=l, nbody 
hmin=hh ( i ) 

sumW ( i ) =mass (i)*W(i,i)/den(i) 
DO jcursor=l, neighb(i) 

j=neighblist (i, jcursor) 

hmin= (hh (i) +hh ( j) ) 12 . 

sumW (i) =sumW (i) tmass (j)*W(i,j)/den(j) 
END DO 
END DO 

! thirdly normalized density 
DO i=l, nbody 

den (i) =den (i) / sumW (i) 
END DO 
END SUBROUTINE DENSUM 
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SUBROUTINE DIVVELO 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

This subroutine computes velocity divergence for particle i 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE ' par am. inc 1 

REAL vji (dim) , rij (dim) , rijO, gradW 
divvel=0 . 
DO i=l, nbody 

DO jcursor=l, neighb (i) 
j=neighblist (i, jcursor) 
hmin= {hh (i) +hh ( j) ) 12 . 
rijO=0 .0 
DO k=l, dim 

vji { k ) =ve 1 ( j , k ) -ve 1 ( i , k ) 
rij (k) =pos (i, k) -pos ( j, k) 
rij0=rij0+rij (k) *ri j (k) 
END DO 

rij 0=SQRT (rijO) 
vdotdelW=0 . 
DO k=l, dim 

gradW-dW (i,j)*rij(k)/rijO 

vdotdelW=vdotdelW+v ji (k) *gradW 
END DO 

divvel ( i ) =divvel ( i ) +mass ( j ) *vdotdelW 
END DO 

divvel ( i ) =divvel { i ) /den ( i ) 
END DO 
END SUBROUTINE DIVVELO 
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SUBROUTINE PRESOUNDENG 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

This subroutine estimates the pressure, sound speed, and 

temperature of particles 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE 'param. inc ' 

! molecular weight relative to the mass of hydrogen 
xmu=2 . 

! hydrogen mass 
xmh=l . 67e-27 

! Bolt zman constant 
xkb=l . 38e-23 

xKK= (xkb/ (xmu*xmh) ) / (ul/ut) **2 
IF (isorad == 1) THEN 

! polytropic index {adiabatic case) 

gamma=2 . 

pre= {gamma-1 ) *den*u 
temp= (gamma-1 ) *u/ xKK 
sound=SQRT (gamma*xKK*temp) 
ELSEIF (isorad == 2) THEN 

! polytropic index (isothermal case) 
gamma=l . 

pre=gamma*xKK*temp*den 
END IF 

END SUBROUTINE PRESOUNDENG 



//////////////////////////////////////////////////////////////////// 



SUBROUTINE RATES 

! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
! This subroutine computes time -rate of velocity (acceleration) 

! and energy of particles 

! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE 'param. inc ' 

REAL rij (dim) , gradW, vij (dim) , mui j 
DO i=l, nbody 

DO k=l, dim 
acc (i, k) =0 . 

END DO 

udot (i) =0 . 

DO jcursor=l, neighb (i) 
j=neighblist (i, jcursor) 
hmin= (hh (i) +hh ( j) ) /2 . 
denmin= (den (i) +den ( j) ) /2 . 
vsig= (sound (i) + sound ( j ) ) 12 . 
ri j = . 
vi jri j=0 . 
DO k=l, dim 
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vi j (k) =vel ( i , k) -vel ( j , k) 
rij (k) =pos (i,k) -pos (j,k) 
ri jO=ri jO+ri j (k) **2 
vi jri j=vi jri j+vi j (k) *ri j (k) 
END DO 

rij 0=SQRT (ri jO) 

IF(vijrij < 0.0) THEN 

mui j=vi jri j*hmin/ (ri j0**2+ (eta*hmin) **2) 

phii j= { -alphas* vsig*mui j+bet as *mui j * *2 ) /denmin 

ELSE 

phii j=0 . 

ENDIF 

preden=pre (i) /den (i) **2+pre (j) /den (j) **2 
rdotdelW=0 . 
vdotdelW=0 . 
DO k=l, dim 

gradW=dW (i, j)*rij(k)/rijO 

acc (i, k) =acc (i, k) -mass ( j ) * (preden+phii j ) *gradW 
rdotdelW=rdotdelW+ri j (k) *gradW 
vdotdelW-vdotdelW+vi j (k) *gradw 
END DO 

udot ( i) =udot (i)+0.5*mass(j)* (preden+phii j ) *vdotdelW 
END DO 
END DO 
END SUBROUTINE RATES 
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FUNCTION W(i, j) 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

This function evaluates the kernel of particles i and j 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE 'param.inc' 
REAL sig, rij, q, W 
ri j=0 . 
DO k=l, dim 

rij=rij+ (pos (i, k) -pos (j, k) ) **2 
END DO 

ri j=SQRT (rij) 
htest= (hh(i) +hh( j) ) 12 . 
IF(hmin /= hh(i) .AND. hmin /= htest) 
& CALL TERROR (' KERNEL : hmin is inconsistent') 

q=ri j /hmin 
! Gauss kernel 
IF (skf == 1) THEN 

sig=( 1.0/3. 14) ** (dim/2.0) 
IF (q <= 3 . 0) THEN 

W=EXP (-q*q) 
ELSE 

W=0 . 
ENDIF 

W-s ig*W/ hmin ** dim 
! spline-base kernel 
ELSEIF (skf == 2) THEN 

IF(dim==l) sig=2. 0/3.0 
IF(dim==2) sig=10 . 0/ (7 . 0*3 . 14) 
IF (dim ==3) sig=1.0/3.14 
IF (q <= 1 . 0) THEN 

W=l .0-1 . 5*q**2+0 . 75*q**3 
ELSEIF (q <= 2.0) THEN 

W=0.25* (2.0-q) **3 
ELSEIF (q > 2 . ) THEN 

W=0 . 
ENDIF 

W-s ig*W/ hmin * * dim 
! Quintic kernel 
ELSEIF (skf == 3) THEN 

IF(dim==l) sig=l. 0/120.0 
IF(dim==2) sig=7 . 0/ (480 . 0*3. 14) 
IF (dim ==3) sig=1.0/ (120.0*3.14) 
IF (q <= 1 . 0) THEN 

W= (3 . 0-q) **5-6. 0* (2 . 0-q) **5+15 .0* (1 .0-q) **5 
ELSEIF (q <= 2.0) THEN 

W= (3 . 0-q) **5-6. 0* (2 . 0-q) **5 
ELSEIF (q <= 3.0) THEN 

W= (3. 0-q) **5 
ELSEIF (q > 3 . ) THEN 

W=0 . 
ENDIF 

W-s ig*W/ hmin * * dim 
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END IF 
END FUNCTION W 
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FUNCTION dW (i, j) 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

This function evaluates the differential of kernel 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INCLUDE 'param.inc' 
REAL sig, rij, q, dW 
ri j=0 . 
DO k=l, dim 

rij=rij+ (pos (i, k) -pos (j, k) ) **2 
END DO 

ri j=SQRT (rij) 
htest= (hh (i) +hh { j) ) /2.0 
IF(hmin /= hh(i) .AND. hmin /= htest) 
s CALL TERROR (' KERNEL : hmin is inconsistent') 

q-ri j /hmin 
! Gauss kernel 
IF (skf == 1) THEN 

sig=( 1.0/3. 14) ** (dim/2.0) 
IF (q <= 3 . 0) THEN 

dW=- (2 . 0*q*EXP (-q*q) ) /hmin 
ELSE 

dW=0 . 
ENDIF 

dW=sig*dW/ hmin** dim 
! spline-base kernel 
ELSEIF (skf == 2) THEN 

IF (dim ==1) sig=2. 0/3.0 

IF (dim ==2) sig=10 . 0/ (7 . 0*3 . 14) 

IF (dim ==3) sig=1.0/3.14 

IF (q <= 1 . 0) THEN 

dW= (-3.0*q+2.25*q**2) /hmin 
ELSEIF (q > 1.0 .AND. q <= 2 . ) THEN 

dW=-0 . 75* (2 . 0-q) **2/hmin 
ELSEIF (q > 2.0) THEN 

dW=0 . 
ENDIF 

dW-sig*dW/ hmin** dim 
! Quintic kernel 
ELSEIF (skf == 3) THEN 

IF (dim ==1) sig=l . 0/120 . 
IF (dim ==2) sig=7.0/ (480.0*3.14) 
IF (dim ==3) sig=1.0/ (120.0*3.14) 
IF (q <= 1 . 0) THEN 

dW= (-5.0* (3. 0-q) **4+30.0*(2.0-q)**4-7 5.0*(1.0-q)**4) /hmin 
ELSEIF (q > 1.0 .AND. q <= 2 . ) THEN 

dW=(-5.0* (3. 0-q) **4+30.0*(2.0-q)**4) /hmin 
ELSEIF (q > 2.0 .AND. q <= 3 . ) THEN 

W=(-5.0* (3. 0-q) **4) /hmin 
ELSEIF (q > 3.0) THEN 

dW=0 . 
ENDIF 

dW-sig*dW/ hmin** dim 
ENDIF 
END FUNCTION dW 
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SUBROUTINE TERROR (message) 

! This subroutine stops the program if there is any error 

I cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
CHARACTER (*) message 

WRITE(*,*) '*********!? i 7 i 7 i error:!?!?!?!*********' 
WRITE)*,*) message 

WRITE!*,*) ' *********program is terminated*********' 
STOP 

END SUBROUTINE TERROR 
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param. inc 



this file contains common definitions and parameters 
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! dimension and maximum number of particles 
INTEGER dim 

! dim — > number of spatial dimensions (1, 2 , or 3) 
PARAMETER (dim=l) 
INTEGER maxn 

! maxn--> maximum number of SPH particles 
PARAMETER (maxn=5 ) 
INTEGER nbody 

! nbody- -> number of SPH particles 
REAL ul, ut, um 

! ul— > unit of length (m) 

! ut— > unit of time ( s ) 

! um— > unit of mass ( kg) 

! ud— > unit of density (kg/m A 3) 

! =um/ul A 3 

! uv— > unit of velocity (m/s) 
! =ul/ut 

! ub— > unit of magnetic field (tesla) 

! =SQRT (um/ul) /ut 

! GO— > gravitational constant 

! =(6.68e-ll*um/ul) * (ut/ul) **2 
COMMON /main/ nbody, ul, ut, um 
! switches for different scenarios 
INTEGER isorad 

! isorad — > isothermal or adiabat ic shock? 

! =1 : adiabatic 

! =2 : isothermal 

INTEGER skf 

! skf — > smoothing kernel function 

! =1: Gauss kernel {Gingold & Monaghan 1981) 

! =2 : spline-base kernel (Monaghan 1985) 

! -3: Quintic kernel (Morris 1997) 
INTEGER nnssl 

! nnssl--> nearest neighbor search and smoothing length 
! =1: tree with h=hf ac* (mass/ den) * * (1/ dim) 

! =2 : tree with dh/dt= (h/dim) *divvel 

! =3 : tree with fixed neighbors between max and min 

INTEGER dsm 

! dsm — > density summation method 

! =1 : summation model without continuity 

! =2: use continuity equation 
COMMON /senarl/ isorad, skf, nnssl, dsm 
! tolerance and correction parameters 
REAL alphas , betas 

! alphas — > shear viscosity 

! betas--> bulk viscosity 
REAL epsi, eta 

! epsi--> parameter in XSPH correction of velocities 

! eta — > parameter to avoid singularities in viscosity 
PARAMETER (epsi=0 .5, eta=0 .1) 
REAL theta, eps 

! theta— > tolerance parameter in tree structure 

! eps — > tolerance parameter in tree structure 
PARAMETER (theta=0 .25, eps=l . Oe-4) 
REAL cour 

! cour--> Cour ant number in step-time 

PARAMETER (cour=0 .25) 

COMMON /toleran/ alphas, betas 

! tree structure data arrays 

INTEGER nsubc, nquad 

! nsubc- -> number of descendants per cell 

! nquad- -> number of independent quadrupole components 

PARAMETER ( nsubc=2 * *dim, nquad=2*dim-l) 

INTEGER inode 

! inode --> initial number of nodes (bodies + cells) 

PARAMETER ( inode=maxn+nsubc*maxn) 

INTEGER mxcell, node, incell, ncell 

! mxcell — > number of cells in the system { =n s ubc* nbody ) 
! node — > number of nodes (bodies + cells) 
! incell — > index of first cell in arrays (=nbody+l ) 
! ncell— > number of cells currently in use (<=mxcell ) 

INTEGER subp (inode, nsubc) , root, sortind (maxn*nsubc) 
! subp--> descendent of each cell 

! root— > index of cell representing root (=incell) 

! sortind (maxn*nsubc) --> sorted cells in descending size 
REAL mid { inode , dim) , clsize (inode) 

! mid — > geometric center of each cell 

! clsize — > size of each cell 
REAL rcrit2 (inode) , quad (inode, nquad) 

! rcrit2 (incell: node )--> critical distances^ of each cell 
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! quad (incell:node, nquad) --> quad moments of each cell 
COMMON /t reel/ mxcell, node, incell, ncell , sort ind 
COMMON /tree 2/ rcrit2, quad, subp, root, mid, clsize, ndesc 
! neiqhbor search parameters and smoothing length 
INTEGER neighb (maxn) , neighblist (maxn, maxn ) 

! neighb (maxn) -> number of neighbors for each particle 

! neighblist (maxn, maxn) — > list of neighbors 
REAL hh (maxn) , hhdot (maxn) 

! hh (maxn) — > smoothing lengths of SPH particles 

! hhdot (maxn) — > smoothing length rate 
REAL hmin 

! hmin--> mean smoothing length of two neighbor particle 
COMMON /neighborl/ neighb, neighblist, hh, hhdot, hmin 
! states of SPH particles 

REAL pos { inode , dim) , mass ( inode) , den (maxn ) 

! pos (node, dim) --> positions of bodies and cells 

! mass (node) --> mass of bodies and cell 
! den (maxn) — > density at position of each particle 
REAL vel (maxn, dim) , divvel (maxn) 

! vel (maxn, dim) --> velocities of each body 
! divvel (maxn) — > divergence of velocity 

REAL acc (maxn, dim) , dendot (maxn) 

! acc (maxn, dim) --> acceleration of bodies 

! dendot (maxn) — > density rate 
REAL sound (maxn) , pre (maxn) , temp (maxn) 

! sound (maxn) --> sound speed of particles 

! pre (maxn) — > pressure of particles 

! temp (maxn) --> temperature of particles 
REAL u (maxn) , udot (maxn) 

! u (maxn) --> energy of particles 

! udot (maxn) --> energy rate 
COMMON / state 1 / pos, mass , vel, divvel , acc, dendot , u, udot 
COMMON /state2/ den, sound, pre, temp 
! time integration parameters 
REAL tnow, dtmin 

! tnow — > current time 

! dtmin — > minimum time-step 
COMMON /time/ tnow, dtmin 



